Magnetic critical behavior of two-dimensional 
random-bond Potts ferromagnets in confined geometries 
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We present a numerical study of 2D random-bond Potts ferromagnets. The model is studied 
both below and above the critical value Qc = 4 which discriminates between second and first-order 
transitions in the pure system. Two geometries are considered, namely cylinders and square-shaped 
systems, and the critical behavior is investigated through conformal invariance techniques which were 
recently shown to be valid, even in the randomness-induced second-order phase transition regime 
Q > 4. In the cylinder geometry, connectivity transfer matrix calculations provide a simple test to 
find the range of disorder amplitudes which is characteristic of the disordered fixed point. The scaling 
dimensions then follow from the exponential decay of correlations along the strip. Monte Carlo 
simulations of spin systems on the other hand are generally performed on systems of rectangular 
shape on the square lattice, but the data are then perturbed by strong surface effects. The conformal 
mapping of a semi-infinite system inside a square enables us to take into account boundary effects 
explicitly and leads to an accurate determination of the scaling dimensions. The techniques are 
appfied to different values of Q in the range 3-64. 



I. INTRODUCTION 



The presence of impurities can have significative effects 
on the nature of phase transitions. Both from experimen- 
tal and theoretical perspectives, the study of the influ- 
ence of randomness is of great importance. Experimen- 
tal evidences of the effect of random quenched impurities 
in two-dimensional systems were found in order-disorder 
phase transitions of adsorbed atomic layers belonging, in 
the pure case, to the Q — 4-state Potts model univer- 
sality class ||l|,|[. In the presence of disorder, the critical 
exponents are modified. On the other hand, no modiflca- 
tion is found when the pure system belongs to the Ising 
universality class Q. 

The study of disordered systems is a quite active field 
of research in statistical physics, and the resort to large- 
scale Monte Carlo simulations is often helpful |j]. Nu- 
merical investigations of the critical properties of ran- 
dom systems require averages over disorder realizations. 
Standard techniques, like Finite-Size Scaling (hereafter 
referred to as FSS) or temperature dependence of the 
physical quantities were extensively used, and, more re- 
cently, conformal invariance techniques were shown to 
provide accurate results. 

The effect of quenched bond randomness in a system 
which undergoes a second-order phase transition in the 
homogeneous case has been considered first. It is well un- 
derstood since Harris proposed a relevance criterion for 
the case of fluctuating interactions pj . Disorder appears 
to be a relevant perturbation when the specific heat ex- 
ponent a of the pure system is positive. Since in the 
two-dimensional Ising model (IM) a vanishes due to the 



logarithmic Onsager singularity, this model was carefully 
studied in the 1980s Q. The analogous situation when 
the pure system exhibits a first-order transition was less 
well studied, in spite of the early work of Imry and Wor- 
tis who argued that quenched disorder could induce a 
second-order phase transition M. This argument was 
then rigorously proved by Aizenman and Wehr, and Hui 
and Berker PJy]. In two dimensions, even an infinitesi- 
mal amount of quenched impurities changes the transi- 
tion into a continuous one. 

The first intensive Monte Carlo (MC) study of the ef- 
fect of disorder at a first-order phase transition is due 
to Chen, Ferrenberg and Landau. These authors stud- 
ied the Q = 8-state two-dimensional random-bond Potts 
model (RBPM), which, in the pure case, is known to ex- 
hibit a first-order phase transition when Q > 4, the larger 
the value of Q, the sharper the transition ||l3]. Taking 
advantage of duality, they performed a finite-size scal- 
ing study at the critical point of a self-dual disordered 
system ||ll|,|l^ and definitively showed that the transi- 
tion becomes of second order in the presence of bond 
randomness. Their results, together with other related 
works [|l3|-|l6|, suggested that any two-dimensional ran- 
dom system should belong to the 2D pure IM universality 
class. These results were also coherent with real experi- 
ments |l|. 

In recent p apers, Cardy and Jacobsen used a differ- 
ent approach |l7[[lq ], based on the connectivity transfer 
matrix (TM) formalism of Blote and Nightingale p9[ . 
They studied random-bond Potts models for different val- 
ues of Q and with a bimodal probability distribution of 
coupling strengths. Their estimations of the critical ex- 



ponents lead to a continuous variation of [S/v with Q. 
This result is in accordance with previous theoretical 
calculations and MC simulations when Q < 4 [pO 21 1. 
In the randomness-induced second-order phase transi- 
tion regime Q > 4, P/i^ is quite different from the Ising 
value of i and particularly in sharp disagreement with 
the Monte Carlo results of Ref. (l|] for Q = 8. Since 
then, Monte Carlo simulations were performed by differ- 
ent groups at Q = 8 ]2^-[24|. The choice of the value 
Q = 8 was motivated by the value of the correlation 
length in the pure case (^ = 23.87 in lattice spacing 
units) 1^. MC simulations which enable to discriminate 
between a first-order regime and a second-order transi- 
tion can indeed be performed easily with systems of larger 
sizes. These studies led to partially conflicting results 
given in Table H, but they eventually found an explana- 
tion in terms of a crossover behavior in a recent work 
of Picco [E4|. While theoretical calculations are gener- 
ally managed in the weak disorder regime (perturbation 
expansion around the homogeneous system fixed point), 
the range of disorder amplitude must be chosen carefully 
in numerical studies, since the random fixed point (FP) 
can be perturbed by crossover effects due to the pure 
and/or the percolation unstable fixed points. The disor- 
dered FP properties are thus more easily observed with 
strong randomness. A disorder amplitude r, given by the 
ratio of the two types of couplings (distributed accord- 
ing to a binary distribution), in the range 8-20 appears 
to be adapted to a numerical analysis and gives a good 
estimate of the disordered fixed point exponents |24|j2q| 
as already observed in the 2D random-bond Ising model 
(RBIM) [|^,|§. 

TABLE I. Bulk magnetic scaling index obtained by differ- 
ent groups in the 8-state Potts model. 



Authors 



0/iy 



Technique 



Chen et al, Ref. [|]J 
Cardy and Jacobsen, Ref. 
Chatelain and Berche, Ref. 
Picco, Ref. H 
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0.118(2) 


MC 


2 


0.142(4) 


TM 


10 


0.153(3) 


MC 


10 


0.153(1) 


MC 



The surface properties of dilute or random-bond mag- 
netic systems were paid less attention. The whole set of 
bulk and surface critical exponents of a given system is 
determined by the anomalous dimensions of the relevant 
scaling fields which enter the homogeneity assumption 
of the singular free energies |g^. The (1,1) surface of 
the disordered Ising model on a square lattice has only 
recently been investigated through MC simulations by 
Selke et al. pQ| . The critical exponent (3i of the boundary 
magnetization was found to be equal within error bars to 
its value in the pure 2D IM. The surface properties of the 
8-state RBPM were also computed in Ref. 1 2^] . 

In this paper, we are interested in the bulk critical 
behavior of disordered Potts ferromagnets, and in the 
evolution of their properties as the number of states Q 



increases. The Hamiltonian of the model is 

~pn = J2K^A,,<r, (1) 

where the spins can take Q different values and the cou- 
pling strengths between nearest neighbor spins are taken 
from a binary probability distribution 

ViK,,) = pSiK,, ~ rK) -K (1 - p)5{K,j - K) (2) 

with p ~ 1/2, which guarantees the self-duality relation 

(e'-^= - l)(e^= - 1) = Q. (3) 

The value r = 1 corresponds to the pure model and 
r ^ oo to the percolation limit. 

In the present work, following previous studies, we use 
the powerful methods of conformal invariance. Talapov 
et al. studied numerically the critical-point correlation 
functions in the 2D RBIM on the torus [Q and took 
into account the finite-size effects through a convenient 
conformal rescaling ^,^. In the cylinder geometry, 
conformal invariance methods have also been success- 
fully applied. In the two-dimensional RBIM, random- 
ness being a marginally irrelevant perturbation, many 
results have been obtained via these techniques: Confor- 
mal anomaly, correlation decay, gap-exponent relation for 
long strips |34|-|3q] . At randomness-induced second-order 
phase transitions, conformal techniques have also been 
used already ||l3,|l8|,0 and numerical evidences for the 
validity of the conformal covariance assumption for cor- 
relation functions and density profiles were recently re- 
ported Q . It is well known that in disordered spin sys- 
tems, the strong fiuctuations of couplings from sample 
to sample require careful averaging procedures |39-H|. 
For that reason, the study of the probability distribu- 
tions must be performed in order to guarantee that the 
average quantities, which should obey the conformal co- 
variance assumption, are correctly obtained numerically. 
A comparison between grand canonical disorder (GCD) 
and canonical disorder (CD) will also be given. 

The plan of the paper is the following: In Sec ||, we 
present the results of connectivity transfer matrix calcu- 
lations on strips with periodic boundary conditions for 
different values of Q. The order parameter correlation 
function, after disorder average, leads to estimates of the 
magnetic scaling index for different strip sizes. From 
our knowledge in the case Q = 8 [^, it appears that 
these computations are suitable for the determination of 
a convenient disorder amplitude in order to reach the dis- 
ordered FP. At large disorder amplitudes (r ~ 10), the 
behavior of the effective central charge can indeed dis- 
criminate between random and percolation fixed points. 
In Sec. [II, we report Monte Carlo simulations in a square 
geometry with the above-mentioned disorder amplitude. 
The magnetization correlation function and density pro- 
file give access to refined values for the corresponding 
exponents. A discussion of the results is given in Sec- 
tion IV. Attention is paid to take into account the differ- 



ent sources of error for the results reported in this work. 



II. CYLINDER GEOMETRY AND DISORDERED 
FIXED POINT 

A. Free energy and central charge 

In the strip geometry, we used the Blote and Nightin- 
gale connectivity transfer matrix method ||l^. In disor- 
dered systems, transfer operators in the time direction 
do not commute and, as a consequence, the free energy 
density is no longer defined by the largest eigenvalue of 
a single TM, but in terms of the leading Lyapunov ex- 
ponent. For a strip of size L with periodic boundary 
conditions, the leading Lyapunov exponent follows from 
the Furstenberg method |42| : 
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(4) 



where T/c is the transfer matrix between columns fc — 1 
and k and \vq) is a suitable unit initial vector. The free 
energy density is thus given by 



[/o(i)]av - -L-lAo(L), 



(5) 



where [. . .]av denotes the average over disorder realiza- 
tions. 

In the following, we considered canonical disorder, a 
situation in which exactly the same numbers of couplings 
K and rK are distributed over the bonds of the whole 
system of length ~ 10^. This choice contributes to re- 
duce sample fluctuations. This is shown in Fig. [l| where 
the stability of the free energy density is compared to the 
standard grand canonical disorder for different runs up 
to m = 10*^ iterations of the TM. 

In Eq. (0), the disorder average is implicitly performed 
through an infinite number of iterations of the trans- 
fer matrix. In our computations, only a finite num- 
ber TO is used, leading to approximate values denoted 
by Aq (L;to) for different runs labelled by an integer 
i — 1, M. The leading Lyapunov exponent and the cor- 
responding eigenvector, | Aq), obtained after to. = 10^ 
iterations of the TM, are then averaged over M = 48 
independent runs. The average free energy density of 
Eq. (J^) is thus replaced, in the calculations, by 



[/0(i)]a 



-L 



'(-^fAi'>(L;10«))^ 



(6) 



The value M = 48 was chosen in order to guarantee a 
stability of the averaged quantities with a relative error 
better than 10~^ for the free energy density and than 
4 X 10~^ for the components of the corresponding eigen- 
vector. The computations are then performed on strips 
of sizes L = 2 to 8. 



O 
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FIG. 1. Free energy density (up to an additive constant 
InQ) vs m, the number of iterations of the TM for a strip 
of size L = 6 (Q = 8, r = 10) with 5 realizations in grand 
canonical (a) and canonical (b) disorder. 

The numerical investigation of critical properties in 
random systems requires the knowledge of the range of 
disorder amplitude (measured here by the ratio r be- 
tween strong and weak couplings) for which the fixed 
point properties is reached. Outside this regime, strong 
crossover effects perturb the data 1 38 . A convenient dis- 



order amplitude r can be obtained from the behavior of 
the effective central charge, which increases when the sys- 
tem approaches the disordered fixed point in non-unitary 
theories as it seems to be the case in the RBPM ]1^,Q . 
The central charge c is defined by the leading size depen- 
dence of the free energy density, and, since the strip sizes 
are quite small, corrections to scaling must be included: 
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The comparison between successive sizes L and L + I al- 
lows us to define a reduced difference which leads to 



[A/i(i)]a 



_ 6[/o(i)]av-[/o(i + 0]. 



(L + iy 

^AX, 

■K 



L- 



where the reduced parameter A is given by 



A 



{L + iy 



L~ 



(8) 



(9) 



In the thermodynamic limit, the central charge c then 
follows from a linear fit as shown in Fig for strips of 
sizes i = 2 to 8 in the case Q — ?>. We restricted our 
study to integer values of r and the data for the effective 
central charge at different disorder amplitudes are given 



in the Table ||. We observe that the value of c is strongly 
depending on the disorder amplitude: It increases from 
the weak disorder limit up to a maximum value and then 
decreases slowly as r increases. 

The central charge at the random fixed point {i.e. 
the maximal value obtained for an optimal disorder am- 
plitude r*{Q)) is shown in Fig. 0. Assuming a linear 
behavior in In Q [p7| which preserves the Ising value 
c{Q = 2) = 1/2 |lj, one gets 

InQ 



c{Q) 



21n2 



(10) 



whilst the percolation limit leads to c(Q) — ^r In Q Q . 
The two behaviors are shown in Fig. 0. The numerical 
data are in good agreement with Eg. JlOJ and are accu- 
rate enough to consider that the random FP has been 
reached at r* (whose values are coherent with those found 
by Jacobsen and Cardy ||l|l: r*(3) = 7, r''{8) = 9 and 
r*(64) — 10). In the following, the scaling properties will 
be studied at the optimal disorder aniplitudes in con- 
tradistinction with previous papers |r^,|l8[ . 
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FIG. 2. Reduced difference between the free energies at 
different sizes [A/i(L)]av for different values of the disorder 
amplitude r {Q — 3). The central charge is given by the in- 
tercept via a linear fit. The parameter I is defined in Eq. (W. 



TABLE II. Extrapolation of the effective central charge c in the thermodynamic limit for the different values of Q and r. 
At each Q, the larger values of c (written in bold face) correspond to the random fixed point regime with an optimal disorder 
amplitude r* . This maximum is not always located at the same r, as shown in the case Q = 3 and 64 for two different runs. 
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FIG. 3. Central charge at the random fixed point as a func- 
tion of the number of states. The full line corresponds to 
c{Q) = lnQ/(21n2) while the dashed line is the percolation 
limit c{Q) = 5\/31nQ/47r. Error bars are smaller than the 
sizes of the symbols. 



[(G'cr(u))]av, our first aim is to show that, in spite of the 
lack of self-averaging, our numerical experiments lead to 
well-defined averages. 
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B. Probability distribution of the correlation 
function 



For a specific disorder realization, the spin-spin corre- 
lation function along the strip 



{G.{u)) - 



Q-l 



(11) 



where (. . .) denotes the thermal average, is given by the 
probability that the spins along some row, at columns j 
and j -\- u, are in the same state (j and j + u measure the 
position in the longitudinal direction of the strip) : 

(Aoig,(ni;r'n.)d,+«iAo) 

(^-.-.+J = — ,,\^L-iJ ,. . — > (12) 



(Aoini 



j+u-1 



Tfc I An 



where | Aq) is the ground state eigenvector and T'j, is 
the transfer matrix in the extended Hilbert space which 
includes the connectivity with the origin site j. The op- 
erator gj identifies the cluster containing Cj , while dj+u 
gives the appropriate weight depending on whether or 
not (Tj+„ is in the same state as aj. The computation is 
performed with a grand canonical disorder. 

An analysis of the correlation function probability dis- 
tribution is needed in order to ensure that self-averaging 
problems do not alter the mean values |4J] . The method- 
ology that we propose is to deduce the critical behavior 
from the decay of the correlation functions using con- 
formal symmetry. Since conformal covariance assump- 
tion is supposed to be satisfied by average quantities^ i.e. 



FIG. 4. Probability distribution of the correlation function 
after 63436 realizations of disorder for a strip of size L — 6 
{Q — 8, r — 10). The vertical dotted line shows the aver- 
age value [{Go-(30))]av while the long-dashed line shows the 
typical value e[<'"'^-(^°^>l--. 



The probability distribution of the correlation func- 
tion, as shown in Fig. 0, enables us to determine the 
most probable (or typical) value G'™p(m) and the aver- 
age correlation function [{Ga-{u))]av, as well as the av- 
eraged logarithm, [ln(G'^(u))]av at any value of the dis- 
tance u. Compatible behaviors are found for G^^{u) and 



It is a confirmation of the essentially log- 



3[ln(G„(«))]a 

normal character of the probability distribution [M , as 
argued by Cardy and Jacobsen [0. It is thus necessary 
to perform averages over larger numbers of samples for 
[(Gcr(u))]av than for [ln(G^(u))]av to get the same rela- 
tive errors. 



Following Cardy and Jacobsen, since the moments 
of the logarithm of the correlation function are self- 
averaging, a cumulant expansion can then be performed 
to reconstruct [(Gcr(u))]av and to compare to the values 
obtained by averaging directly over the samples. 
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FIG. 5. Average correlation function, most probable (or 
typical) value and sum up to the 4th order of the cumulant 
expansion obtained from 63436 realizations of disorder for a 
strip of size L = & [Q = 8, r = 10). 

The results in Fig. || (for Q = 8), strengthen the credi- 
bility of the direct average, and also clearly show that the 
cumulant expansion up to fourth order still strongly fluc- 
tuates at large distances compared to \{Ga{u))]av In the 
following we will thus favour the direct averaging process, 
using a large number of disorder realizations. 



where u = Re {w2 — wi). The scaling index x^ can thus 
be deduced from a linear fit in a semilog plot. 

For each strip size [L ~ 2 — 8), we realized 80 x 10^ 
disorder configurations. It allowed us to define mean 
values and error bars for the correlation functions at 
any point in the range u = \ — 100, taking into ac- 
count the standard deviation over the samples. The non 
self-averaging behavior of the correlation functions in- 
duces large variances (The reduced variance Rx{L) = 
([X^]av — [X]av)/[-^]av ^ocs uot bchavc as a power law, 
but evolves towards a constant value when the strip 
size increases, e.g. Rq^(2Q){L) -^ 1.50, as already ob- 
served for several quantities by Wiseman and Domany in 
Refs. pO| , [4l|| . ) . The exponents follow from an exponen- 
tial fit in the range w, > 5 and [(Go-(M))]av > e, where 
the cutoff e is introduced in order to avoid tiny numbers 
whose values are lower than the fluctuations. The error 
bars given for the exponents take into account the un- 
certainties of data for the correlation functions |^^ . The 
resulting values for each strip size are plotted against 
L~^ which allows an extrapolation in the thermodynamic 
limit. This is shown in Fig. ^ in the case Q — ^. This 
figure provides a confirmation of the effect of a too weak 
disorder: Strong crossover effects take place which lead 
to a wrong determination of the critical behavior with 
the strip sizes used here. On the other hand, at the opti- 
mal value r*(Q), the exponent converges in the L — > oo 
limit towards a well defined final estimate. 



C. Bulk magnetic scaling dimension 

We will now use the results that follow from the as- 
sumption of conformal covariance of the average correla- 
tion functions. In the infinite complex plane z = x + iy 
(denoted by the index cxd) the correlation function ex- 
hibits the usual algebraic decay at the critical point 



[(G,(i?))]av = [('T(zi)a(z2))oo]av 

— 9 ^ 

= const X R ^^^, 



(13) 



where R =\ zi — Z2 \ and x'^ — (3/v is the bulk magnetic 
scaling dimension. Under a conformal mapping w(z), 
the correlation functions of a conformally invariant 2D- 
system transforms into the new geometry according to 

G.{W^,W2) =1 W'{Z^) r^' I W'{Z2) r"' G,(Z1, Z2). (14) 

The logarithmic tranformation w = -k-lnz is known to 
map the z plane onto an infinite strip (denoted by the 
index st) w = u + iu of width L with periodic bound- 
ary conditions in the transverse direction. Applying 
Eq. (nj) in the random system where [(Gcr(z«i, W2))]iw = 
[(cr(?x;i)cr(ii;2))st]av corresponds to the strip geometry, one 
gets the usual exponential decay along the strip 



[{GAu))]. 



const X exp 



27r 



xtu 



(15) 




FIG. 6. Magnetic scaling index deduced from the algebraic 
decay of the average correlation function along the strip of 
size L as a function of L~^ and extrapolation in the thermo- 
dynamic limit (Q = 8, L = 2 — 9 for r = 2 and 20, from 
Ref. [38] and L = 2 — 8 for r = 10, this work, where the data 
analysis is more refined (see Appendix A) leading to error 
bars 10 times smaller). 



The convergence of effective scaling dimensions at dif- 
ferent strip sizes, obtained with a cutoff value in the range 
e — lO^"' — 10~^ and r = r*{Q), is shown in Fig. |7| for 
different values of Q. The ex trap olation in the thermody- 
namic limit is given in Table 



III 



The details of the fitting 
procedure and of the evaluation of errors is presented in 
Appendix ^. 
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FIG. 7. Magnetic scaling index deduced from the algebraic 
decay of the average correlation function along the strip of 
size L as a function of L~^ and extrapolation in the thermo- 
dynamic limit for different Q-values (L = 2 — 8). 



TABLE in. Bulk magnetic scaling index (after extrapo- 
lation in the thermodynamic limit) obtained from the decay 
of the correlation function along the strip (cutoff parameter 
e = 10"* - 10""). 



Q 


r 
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Ax^ 


3 


5 


0.1321 


3 X IQ-'* 


4 


7 


0.1385 


3 X IQ-* 


5 


7 


0.1423 


3 X IQ-* 


6 


8 


0.1456 


3 X 10-* 


8 


10 


0.1505 


3 X 10~* 


15 


10 


0.1572 


3 X 10-* 


64 


12 


0.1669 


3 X 10"'' 
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FIG. 8. a) Monte Carlo simulations of the 2d RBPM inside 
a square of 101 x 101 lattice sites (10'^ MCS/spin, Swend- 
sen- Wang cluster algorithm) . The figure shows the correlation 
function between a point close to the surface {C,i — i) and all 
other points ( in the square. The notations are specified in 
b). 



III. SQUARE GEOMETRY AND CRITICAL 
BEHAVIOR 



A. Conformal rescaling of boundary effects 

Monte Carlo simulations of two-dimensional spin sys- 
tems are generally performed on systems of square shape. 
In the following, we consider such a system of size Nx N, 
and call u and v the corresponding directions (Fig. H). 



The order parameter correlation function between a 
point close to the surface, and a point in the bulk of the 
system should, in principle, lead to both surface and bulk 
critical exponents, possibly to structure constants ||46| . 
Practically, FSS techniques are not of great help for the 
accurate determination of critical exponents, since 

i) strong surface effects (shape effects) occur which 
modify the large distance power-law behavior, i.e. the 
scaling regime. 



ii) the universal scaling function entering the correla- 
tion function is likely to display a crossover before its 
asymptotic regime is reached (system-dependent effect). 

One can proceed as follows; Systems of increasing sizes 
are successively considered, and the correlations are com- 
puted along u- (parallel to a square edge considered as 
the free surface) and w-axis (perpendicular to this edge). 
The order parameter correlation function for example is 
supposed to obey a scaling form which reproduces the ex- 
pected power-law behavior in the thermodynamic limit: 



G7{^) = 



G-(^ 



1 



,x%+x 



1 J± 



(^). 



,1 J\\ 



(^)' 



(16) 



(17) 



where x^ and x^ are the bulk and surface order parame- 
ter scaling dimensions, respectively. The scaling func- 
tions have to satisfy asymptotic expansions including 
corrections to scaling due to the limitations mentioned 
above, e.g. f^ (f ) -^ 1 + const x (^)^ -f- . . . in the 
boundary region u — > iV. 

Equations dlq) and (O) are not very useful for the de- 
termination of critical exponents, since the scaling regime 
V —^ N is perturbed by the correction terms which have a 
large amplitude, resulting from the significance of finite- 
size corrections. Nevertheless, conformal invariance sup- 
plies an easy way to take into account explicitly shape 
effects in two-dimensional systems, and thus provides 
a refined procedure for the determination of the expo- 
nents. In pure systems, density profiles, correlations and 
local properties have been investigated in various geome- 
tries (surfaces p^H49|, c orners |50-5|], strips |5^-|5^ ] or 
parabolic shapes |l56|-|6l[, for a review, see Ref. |]62|), as 
well as the moments of the magnetization |63|| and struc- 
ture factors Q have been calculated in square systems. 

In the following, we shall consider a square system with 
free or fixed boundary conditions on all the edges. Us- 
ing conformal invariance techniques pq] , the Schwarz- 
Christoffel mapping enables us to calculate the surface- 
bulk correlation fonction inside the square. The mapping 
of the complex half-plane z = x + iy, Im z > 0, inside a 
square C = u + iv, -N/2 < Re C < N/2, <lm ( < N, 
is realized by the conformal transformation |66[| 



dz 



C 



^(l-z2)(l-fc2z2) 



(18) 



Since ( = N/2 and ^ = N/2 + iN are mapped onto z = 1 
and z — 1/fc (0 < fc < 1), respectively, the constant C is 
related to the size of the square 



N/2C = K{k) = K, 
N/C = K(fc') = K', 



(19) 



where k' = \J\ — fc^ and K(fc) is the complete elliptic in- 
tegral of the first kind. The modulus fc also follows from 
these equations. It is given by Bq] 



k 






(p+l/2)' 



-It: 



The complete transformation is finally written 



N 
C= — F(z,fc) 
^ 2K ^ ' ^ 



N 
K 



-F(z,fc), 



K'C /K'C , 



(20) 



(21) 



(22) 



where F(2:, fc) is the elliptic integral of the first kind and 
sn((^, fc) the Jacobian elliptic sine [p7|. 



B. Correlation functions 

The two-point correlation function of a conformally in- 
variant system can now be obtained in the C~geometry 
in terms of its counterpart in the semi-infinite system 
(z— geometry): 



G(Cl,C)^(^(Cl)cT(C))sq 

= ic'(^i)r^^ic'(^) 



((T(zi)a(z))hp, (23) 



where the correlation function in the half-plane (hp) ge- 
ometry is known to take the form [|47| 



G(zi,z) = (cr(zi)cr(z))hp 

= const X (yiy)"^-'i/'(^), 



(24) 



where the dependence on a; = , ^j^^^ of the universal 
scaling function i]) is constrained by the special conformal 
transformation and its asymptotic behavior, tliiiS) '^ w^" , 
in the limit y\ — 0(1), yS> 1, is implied by scaling. 

Equations. ( p4|) and (|2J^), applied in the random situ- 
ation, lead to the correlations between Ci = i, close to a 
side of the square, and any point inside it, as follows: 

[(G,(C))]av = const X {| C!{z) I Im(2(C))}-^^VXc.). (25) 

Taking the logarithm of both sides, the bulk critical ex- 
ponent x\ can thus be deduced from a linear fit along 
•jj — const curves in the square: 



ln[(G,(C))]av =const'-4lnK(C) + ln7^(w), (26) 



with 



^(C)^Im(z(C))|[l-z2(C)][l-fcV(C)] 



1-1/2 



(27) 



We will now discuss the results of MC simulations per- 
formed with the Swendsen-Wang cluster algorithm |6q] 
for systems of size 101 x 101 with canonical disorder. 
The details concerning the choice of the parameters for 
the simulations (number of MC iterations, . . . ) are given 



in Appendix |B[ Average over disorder is performed over 
-^rdm = 3000 samples. All the MC simulations are done 
at the optimal disorder amplitude r*{Q) determined in 
the strip geometry. 

Eq. ( Pq ) is used in Fig. ^ to extract the bulk magne- 
tization scaling dimension at Q = 8. Consistent values 
are obtained for different fixed values of the parameter 
w. Averaging the results at different w's, one obtains 



4(8) = 0.152 ±0.003, 
corresponding to an error of 2%. 
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FIG. 9. Rescaled correlation function along six u = const 
curves in the square (shown in the upper part). These curves 
are approximated by linear expansions in the neighborhood 
of the discrete lattice sites, which explains the variations on 
the sizes of error bars (Q = 8, r = 10). 

One should nevertheless mention that the uncertainty 
on this result is underestimated, since neither the fluctu- 
ations due to randomness, nor the influence of a variation 
of r around the optimal value has been taken into account 
explicitly. This is intentional, since such studies would 
require intensive computational efforts and would be less 
accurate that the next method to be presented. 



C. Density profiles 

Owing to the unknown scaling function ip{uj), the de- 
termination of the bulk critical exponent from the behav- 
ior of the correlation function is not extremely accurate. 
Furthemore, since a few points are used for the flts along 
u = const, curves, this introduces a poor statistics. It can 
nevertheless be improved if one considers the magnetiza- 
tion profile inside a square with fixed boundary condi- 
tions. Since it is a one-point function, its decay from the 
distance to the surface in the semi-infinite geometry is 
fixed, up to a constant prefactor 



[(0'(2;))hp]av -- y~ 



(29) 



The local order parameter is defined, according to 
Ref. ||6^, as the probability for the spin at site C in the 
square, to belong to the majority orientation (Fig. 00). 
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FIG. 10. Density profiles inside the square averaged for 
3000 disorder realizations (Q = 8, r = 10). 

The Schwarz-Christoffel mapping leads to the following 
expression for the average profile in the square geometry: 



[(c^(C))]av = const X 



v/|i-z2(C)| -11-^^^(01 

Im(z(C)) 



(30) 



This expression, of the form [(cr(C))]av = [/(^)/2/]^"! 
holds for any point inside the square. It allows an accu- 
rate determination of the critical exponent, since the N"^ 
lattice points enter the power-law fit (Fig. O). Although 
this technique is more precise than the previous one, one 
has to take care to different sources of error. It is indeed 
again necessary to consider the influence of the number 
of disorder configurations which are used to get the av- 
erage magnetization, as well as the effect of a variation 
of the disorder amplitude around the optimal value. We 
performed A^rdm = 5000 realizations of disorder in five 
independent runs (see Appendix H), and computed the 
magnetic exponent for each run. Averaging the results, 
it yields the values given in Table IV. 
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FIG. 11. Rescaled magnetization profile inside the square 
for 5000 disorder realizations (Q = 3, 8 and 64 from top to 
bottom). The power law fits are over 100^ data points. 



TABLE IV. Bulk magnetic scaling index obtained from the 
magnetization profile inside the square (5000 realizations of 
disorder) . 

Q r x^ Axc, 

3 X 10"^ 

4 X IQ-^ 

4 X 10"^ 

5 X 10~^ 

5 X 10"^ 

6 X 10~^ 
6 X 10^^ 



3 


5 


0.13357 


4 


7 


0.13815 


5 


7 


0.14302 


6 


8 


0.14621 


8 


10 


0.15031 


15 


10 


0.15984 


64 


12 


0.17299 



the case of the bulk, but the estimation x^(8) ~ 0.47(3) 
is in agreement with the value that we obtained previ- 
ously by FSS techniques in Ref. [|2|. It also agrees with 
the TM results which give xl{8) ~ 0.48(2) for L = 7 and 
xl{8) ~ 0.50(2) for L = 8. 
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FIG. 12. Large distance behavior of the universal scaling 
function {Q = 8, r — 10) leading to the surface scaling index. 
The fit has been shifted for clarity. 

If the leading singularity (xl.) is found to be the same 
using the two techniques, we note that the corrections to 
scaling are very different, as it appears in the deviation 
between the curves as uj increases. This can be the result 
of the ensemble average procedure which is not identical 
in the two approaches (grand canonical for the TM tech- 
nique and canonical disorder for the MC simulations). 
The same type of sensitivity to the ensemble average was 
reported recently by Wiseman and Domany Mw. 



D. Boundary critical behavior 

The surface scaling dimension can be obtained once 
the bulk exponent is known. From standard scaling, the 
asymptotic behavior of the two-point correlation func- 
tion, when j/i = 0(1), y S> 1 is expected to involve both 
bulk and surface dimensions: 



[(G.(2/-yi))]av~y-(^-+^^). 



(31) 



A power law behavior thus follows for the universal 
scaling function defined in Eq. (24): 



V'(c^)^[(G.(C))]avX{|C'(z)|Im(z(C))F° 

~ LU^", UJ ^ 0. 



(32) 



A log-log plot of Eq. ( |32D is shown on Fig. ^, where 
the TM results are also presented for comparison. The 
result for the surface scaling index is less accurate than in 



IV. CONCLUSION 

In this paper, we have investigated the magnetic criti- 
cal properties of disordered Q— state Potts ferromagnets 
for a wide range of Q— values. These models lead to 
second-order phase transitions which are particularly in- 
teresting, since they belong to new universality classes. 
The accurate determination of critical indices is a pre- 
liminary step towards a deeper understanding of these 
universality classes. Although universality is expected 
with respect to the disorder amplitude r, previous works 
on finite systems have shown that the numerical results 
are very sensitive to the choice of this disorder amplitude. 
This sensitivity is attributed to crossover effects due to 
the pure model (r — > 1) and percolation (r — > oo) un- 
stable fixed points. The behavior of the effective central 
charge as a function of r can fortunately be exploited 
to locate the optimal regime of disorder. One should 
mention that in our previous studies, this extreme sen- 
sitivity of the numerical estimates of critical exponents 
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was not well understood, resulting in an underestimation 
of uncertainties. We tried to present here a careful! anal- 
ysis leading to reliable error bars. This uncertainty is 
mainly due to the non self-averaging behavior of corre- 
lation functions. In the strip geometry, the number of 
samples being already important, better estimates would 
not be easy to obtain, whilst in the MC simulations, im- 
provements could be supplied by increasing the number 
of realizations of disorder. 

The conformal mapping inside the square seems very 
efficient compared to standard FSS studies, one lattice 
size being needed only. The accuracy is furthermore sub- 
stantially improved, since 

i) the finite-size corrections are essentially included in 
the conformal mapping, 

ii) all the lattice points enter the fit of the density pro- 
files. 

A summary of our results, compared to other inde- 
pendent determinations of the magnetic scaling index, is 
given in Table M and the dependence on Q is shown in 
Fig. O. The pure model value for Q < 4 is shown for 



comparison [Q. Both FSS and conformal invariance re- 
sults are presented. The two techniques used in this work 
are in agreement with each other, as well as with previous 
studies at the same disorder amplitude, at least as long 
as the number of states Q is not too large. When the 
ratio r is very different, disagreement with other studies 
which are likely due to crossover effects occurs. On the 
other hand, when the number of states is large, Q > 15, 
there appears discrepancies between the two techniques 
used here. Whilst the second method (square geometry, 
5000 realizations) seems to be the most accurate, we are 
more confident in the first one (strip geometry, 80000 re- 
alizations): If the number of disorder realisations is too 
small, the average behaviour will indeed give an exponent 
closer to the typical one, and thus too large. MC simula- 
tions are furthermore known to be less efficient when Q 
increases, since the autocorrelation time increases also, 
requiring larger numbers of thermalisation iterations. 

We also note that the leading singularity of the mag- 
netization does not depend, up to the precision of our 
results, on the type of disorder considered (GCD or CD). 



TABLE V. Extrapolation of the bulk magnetic scaling dimension x^ in the thermodynamic limit for the different values of Q. 
The first two columns recall previous FSS results obtained by MC simulations (in which the accuracy had been overestimated, 
since the influence of the disorder amplitude was not well understood, at least in which concerns our own studies). The data in 
the four remaining columns were deduced from conformal invariance. The quantity that was studied is indicated in the table 
as well as the geometry and the numerical technique. The results presented in this work are written in bold face. The table 
notes recall the parameters used for each result, especially the values of disorder amplitude which are known to have strong 
influence on the exponent. 
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0.15984(6)'' 
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='MC simulations (Wolff algorithm, ~ lO'' samples, r = 10, GCD) from Ref. |2ll. 

''TM calculations (L = 1 - 7, 10^ samples, r = 2, GCD) from Refs. ^,^. 

'=TM calculations (L = 2 — 8, 80 x 10'^ samples, the values of disorder amplitude for Q = 3, 4, 5, 6, 8, 15 and 64 are r = 5, 7, 

7, 8, 10, 10 and 12, respectively, GCD), this work. 

"^MC simulations (Swendsen-Wang algorithm, N = 101, 5 x 10'^ samples, the values of disorder amplitude for Q = 3, 4, 5, 6, 8, 

15 and 64 are r = 5, 7, 7, 8, 10, 10 and 12, respectively, CD), this work. 

"MC simulations (cluster algorithm, N = 256, ~ 500 samples, r = 10, GCD) from Ref. ^. 

''MC simulations (Wolff algorithm) M. Picco, Ref. [52] cited in Ref. ||l||. 

^MC simulations (Swendsen-Wang algorithm, N < 100, ~ 30 samples, r = 2, restricted CD) from Ref. |12|. 



''MC simulations (Wolff algorithm, A^ < 100 and 500, ~ 10^ samples, r = 10, GCD) from Ref. 14 1 

500 samples, r = 10, CD) from Ref. [gS] 
10, GCD) from Ref. ^. 
101, 3 X 10^ samples, r = 10, CD) from Ref E 



'MC simulations (Swendsen-Wang algorithm, N < 100, 
JTM calculations (L = 2 - 9, 40 x 10^ samples, r 
""MC simulations (Swendsen-Wang algorithm, N 
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FIG. 13. Q-dependence of the bulk magnetic scaling di- 
mension in the RBPM compared to the pure model value for 
Q<4. 

In the case Q = 3 which was already considered by dif- 
ferent authors, there exists a perturbative result (renor- 
malization group approach for the perturbative series 
around the pure model conformal field theory): 

xl = ^ + 0.00132 ~ 0.13465 (33) 

15 

This result was confirmed numerically by Picco | pl[ and 
Cardy and Jacobsen |17| . In this work, we obtain a value 



which is slightly too small. We nevertheless note that the 
two values, at r = 5 and r = 6 are in perfect agreement 
(see Appendix |b|) . 

We eventually mention a recent work of Olson and 
Young [Q who performed a MC study of the multiscal- 
ing properties of the correlation functions for different 
values of Q. They used a different self-dual probability 
distribution of the couplings, and obtained slightly dif- 



ferent results (e.g. 



0.161(3) at Q = 
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APPENDIX A: EVALUATION OF ERRORS IN THE TRANSFER MATRIX CALCULATIONS 



In spite of the large number of disorder realizations, 
the correlation functions along the strip display an im- 
portant dispersion but the resulting values for the critical 
exponents are extremely accurate. In order to obtain a 
correct estimation of the errors on the magnetic scaling 
index, we studied the influence of the cutoff parameter e. 
For e ~ 10~^, a few points are taken into account only 
and the short distance behavior of the correlation func- 



tion is observed. On the other hand, with e ~ 10~^, all 
the data points in the range u = 5 — 100 are taken into 
account in the fit, giving a greater weight to the long- 
distance behavior. Clearly, one has to find a compromise 
between the two approaches. Fortunately a variation of 
the cutoff parameter does not affect the value of the ex- 
trapolated exponent which remains very stable as shown 
in Table W\ 
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TABLE VI. Bulk magnetic scaling index (after extrapolation in the thermodynamic limit) obtained from the decay of the 
correlation function along the strip with different values of the cutoff e. 



Effective exponent at Q = 3, r : 



0.125 xlO" 

0.13207 

5.3x10"* 



0.312 xlO" 

0.13209 

3.8x10"'' 



0.781 X 
0.13209 
3.3x10" 



10 



0.195 xlO" 

0.13208 

3.1x10"** 



0.488 XlO" 

0.13209 

3.1x10"'' 



0.122 XlO" 

0.13209 

3.0x10"** 



0.305 XlO" 

0.13209 

3.0x10"" 



0.763 XlO" 

0.13209 

3.0x10"" 



Effective exponent at Q = 8, r = 



10 



^ 0.763 XlO"" 
0.15054 
2.6x10"" 



Ax 



0.125 XlO" 

0.15014 

4.7x10"** 



0.312 xlO" 

0.15032 

3.4x10"" 



0.781 X 
0.15047 
2.9x10" 



10 



0.195 XlO" 

0.15050 

2.7x10"" 



0.488 XlO" 

0.15050 

2.7x10"" 



0.122 XlO" 

0.15053 

2.6x10"" 



0.305 XlO" 

0.15054 

2.6x10"" 



Effective exponent at Q = 8 
"J 



11 



^ 0.763 XlO"" 
0.15078 
2.6x10"" 






0.125 xlO" 

0.15040 

4.7x10"" 



0.312 xlO" 

0.15056 

3.3x10"" 



0.781 X 
0.15072 
2.9x10" 



10 



0.195 XlO" 

0.15071 

2.7x10"" 



0.488 xlO" 

0.15074 

2.6x10"" 



0.122 xlO" 

0.15077 

2.6x10"" 



0.305 XlO 

0.15078 

2.6x10"" 



Effective exponent at Q 
7 



64, 



12 



"^ 0.763 XlO"" 
0.1671 
3x10"" 



Ax 



0.125 XlO" 
0.1663 

6x10"" 



0.312 XlO" 

0.1663 

4x10"" 



0.781 X 

0.1667 

3x10"* 



10 



0.195 XlO" 

0.1668 

3x10"" 



0.488 XlO" 

0.1668 

3x10"* 



0.122 XlO" 

0.1669 

3x10"* 



0.305 XlO 

0.1670 

3x10"" 



Another contribution to the error should come from 
the choice of the disorder amplitude. To study this ef- 
fect, we considered a variation of r close to the optimal 
value. It leads to a result which is inside the error bars of 
the previous one, as shown in the case Q = 8 in Table [V]|. 
The uncertainty in the range e = 10~* — 10"^ is of the 
same order of magnitude than the fluctuations between 
the data obtained with different values of e and r, so we 
eventually consider as a definitive result the fit with this 
cutoff value. 



APPENDIX B: DETAILS OF THE MONTE 
CARLO SIMULATIONS 

In random systems, in addition to the usual MC error, 
the random-bond fluctuations introduce another source 
of statistical error. For any physical quantity X, the total 
error is given by 



(SX)' 



rdin 



4(H-2tx) 



(Bl) 



Nidm NrdmNuC 

where the first term is due to the disorder fluctuations, 
whilst the second one describes the fluctuations during 
the MC iterations. This latter term corresponds to the 
standard deviation of independent random variables, cor- 
rected by the autocorrelation time to take into account 
the correlations between the successive data. In these ex- 
pressions, A'^MC is the number of MC iterations, measured 
in MC steps (MCS), realized for the measurements of the 
physical quantities for each disorder realization, Nj-^m is 
the number of disorder realizations and tx is the auto- 
correlation time for the quantity X (the definition of tx 



sometimes absorbs the factor 1 describing uncorrelated 
variables). The variances ctt and (Trdm respectively mea- 
sure the deviation due to thermal fluctuations for a given 
sample and the deviation from the exact value within the 
ensemble of disorder configurations. 

Both variances are of the same order of magnitude. 
The leading source of error thus comes from the disorder 
average and a large number of samples is needed in order 
to get accurate results. In our simulations we used the 
Swendsen-Wang cluster algorithm JGgl for systems of size 
101 x 101. The autocorrelation time for the total magne- 
tization is Tcr ~ 35 MCS. The prehminary 5000 MCS have 
been discarded for thermalisation (better that 10^ x Tu), 
and A'^MC — W^ MCS were done to compute the physi- 
cal quantities. Average over disorder is performed over 
-^rdm — 5000 samples. From preliminary runs over 1000 
samples, we deduced the standard deviations a-^^m — ^■^'^ 
and cr^Q ~ 0.13. The order of magnitude of the two con- 
tributions to the error is thus 



5a 



MC 



Sa, 



rdn 




4(1±2Z^ , 6 X 10"^ 



^li^c^ 1.36 X 10"^ 

Nidni 



(B2) 



for a point in the middle of the square. Due to the flxed 
boundary conditions, close to the edges of the square the 
fluctuations are reduced. These values conflrm the signif- 
icance of the disorder contribution (tt^tF- — 0.08% and 

r ;^W'" ~ 2%). The values of the exponent x'^ for different 



values of Q and r are given in Table VIl 
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TABLE VII. Bulk magnetic scaling index obtained from the profile of the order parameter inside a square with fixed 
boundary conditions for five independent runs (1000 configurations of disorder for each run). The final result obtained with 
5000 configurations of disorder is given in the column called average. 









Exponent 


at Q 


= 3, r = 5 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


xi 


0.13260 


0.13384 


0.13405 




0.13418 


0.13323 


0.13357 


^xi 


7x10^5 


7x10-5 


7x10-5 




7x10-5 


7x10-5 


3x10-5 








Exponent 


at Q 


= 3, r = 6 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


xi 


0.13450 


0.13262 


0.13307 




0.13378 


0.13333 


0.13345 


^x\ 


7x10^^^ 


7x10-5 


7x10-5 




7x10-5 


7x10-5 


3x10-5 








Exponent 


at Q 


= 4, 7- = 7 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


4 


0.13886 


0.13835 


0.13798 




0.13703 


0.13858 


0.13815 


Ax^ 


8xl0~5 


8x10-5 


8x10-5 




8x10-5 


8x10-5 


4x10-5 








Exponent 


at Q 


= 4, r = 8 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


4 


0.13799 


0.13794 


0.13753 




0.13849 


0.13798 


0.1379 


Aa;* 


9x10^5 


9x10-5 


9x10-5 




9x10-5 


9x10-5 


4x10-5 








Exponent 


at Q -- 


= 8, r = 10 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


xi 


0.1508 


0.1515 


0.1501 




0.1501 


0.1492 


0.15031 


^x\ 


IxlO""* 


IxlQ-"* 


1x10-'' 




IxlQ-"* 


1x10-'' 


5x10-5 








Exponent 


at Q -- 


= 8, r = 20 








run 1 


run 2 


run 3 








average 


xl 


0.14527 


0.14506 


0.14505 








0.14513 


Ax^ 


6x10"^ 


6 X 10-5 


6 X 10-5 








3 X 10-5 








Exponent at Q = 


= 64, r = 12 








run 1 


run 2 


run 3 




run 4 


run 5 


average 


4 


0.1722 


0.1733 


0.1724 




0.1747 


0.1725 


0.17299 


Ax^ 


1 X 10"^ 


1 X 10-^ 


1 X 10"'' 




1 X 10"* 


1 X 10-" 


6 X 10-5 



' To whom all correspondence should be addressed. Elec- 
tronic address: berche@lps.u-nancy.fr 
The Laboratoire de Physique des Materiaux is Unite 
Mixte de Recherche C.N.R.S. No. 7556. 
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